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Abstract 

Temporal integration of equations possessing continuous symmetries (e.g. systems with translational invariance as- 
sociated with traveling solutions and scale invariance associated with self-similar solutions) in a "co-evolving" frame 
(i.e. a frame which is co-traveling, co-collapsing or co-exploding with the evolving solution) leads to improved ac- 
curacy because of the smaller time derivative in the new spatial frame. The slower time behavior permits the use 
of projective and coarse projective integration with longer projective steps in the computation of the time evolution 
of partial differential equations and multiscale systems, respectively. These methods are also demonstrated to be 
effective for systems which only approximately or asymptotically possess continuous symmetries. The ideas of pro- 
jective integration in a co-evolving frame are illustrated on the one- dimensional, translationally invariant Nagumo 
partial differential equation (PDE) . A corresponding kinetic Monte Carlo model, motivated from the Nagumo kinet- 
ics, is used to illustrate the coarse-grained method. A simple, one-dimensional diffusion problem is used to illustrate 
the scale invariant case. The efficiency of projective integration in the co-evolving frame for both the macroscopic 
diffusion PDE and for a random-walker particle based model is again demonstrated. 

Key words: projective integration, coarse projective integration, continuous symmetry, multiscale computation, dynamic 
renormalization 



1. Introduction 



Projective and coarse projective integration have been recently proposed as effective methods for the 
computation of long time behavior in complex multiscale problems ^J, ^1 UtI . The main idea is to 
use short bursts of appropriately initialized simulations to estimate the time derivative of the quantities of 
interest and then use polynomial extrapolation to jump forward in time |l3lll8j| . When projective integration 
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is applied to deterministic problems (governed by systems of differential equations), one can show that 
it might significantly accelerate the computation of time evolution for systems with large gaps in their 
eigenvalue spectrum [16J. By wrapping the same algorithm around an inner atomistic and/or stochastic 
simulator, one can similarly accelerate coarse-grained computations fl^ IT3 . flij ] . 

Many problems possess additional continuous symmetries 0, flfl l2fl l25l| which can give rise to solutions 
which are traveling, exploding, collapsing or rotating in the domain of interest. In principle, projective 
integration might be applied to such systems as well. However, we can improve the efficiency of the method 
by taking the underlying symmetry into account. The key idea is to perform the projective integration in a 
"co-evolving" frame |32l ]. 

Projective integration in a co-traveling frame is applied to the Nagumo equation [2H Esj |. a well-studied 
system with translational invariance and traveling solutions. Projecting in a dynamically renormalized frame 
is similarly applied to systems characterized by scale invariance. The scale invariant system we study in this 
paper is one-dimensional diffusion. In both applications the existence of continuous symmetries is exploited. 
We apply modified projective integration protocols that are implemented in a dynamically co-evolving frame. 
We demonstrate that this modification improves computational accuracy, allowing for large projective steps. 

We also illustrate the coarse-grained version of projective integration in a co-traveling frame, for a Stochas- 
tic Simulation Algorithm (SSA) [l<| implementation of the Nagumo kinetics. The density of reactant particles 
progressively forms a traveling wave front moving with a constant shape and velocity. In a second application 
we study one-dimensional diffusion simulated by a large ensemble of random walkers. The macroscopic behav- 
ior, described by the cumulative density function (CDF) of the particle positions, features scale invariance. 
Projective integration is appropriately modified to exploit this scale invariant character of the macroscopic 
behavior and increase the accuracy of computations, again allowing for relatively large projective time steps. 

The paper is organized as follows. In Section 2 we briefly summarize projective integration techniques 
and discuss the general ideas of equation-free techniques |23| , a computational framework wrapped around 
microscopic (e.g. kinetic Monte Carlo) simulators. We then present our projective and coarse projective 
integration scheme in a co-evolving frame and its application to problems with translational invariance 
(Section 2.2) as well as scale invariance (Section 2.3). In Section 3 we illustrate the efficiency of the method 
in a co-traveling frame for the Nagumo equation. We also describe coarse projective integration for a kinetic 
Monte Carlo simulation of a reaction-diffusion system based on Nagumo kinetics. In Section 4 we present 
results from the application of projective integration to the scale invariant diffusion system, both for the 
macroscopic diffusion equation and for a random-walker model in one spatial dimension. In Section 5 we 
propose a more general approach that can handle systems evolving in space and scale with an asymptotically 
invariant form, and summarize our work in Section 6. 



2. Projective and Coarse Projective Integration 

Consider a system described by either a suitable macroscopic evolution equation or a stochastic, individual- 
based model, and let M{t) be the macroscopic observable, for which a closed macroscopic evolution equation 
exists. Depending on the problem, M(t) can be a single scalar, a vector or a point in a suitable infinite- 
dimensional Banach space (e.g. a function of physical space). In our illustrative numerical examples the 
macroscopic observable will be a (discretized) field of the density of individuals. 

The parameters of the forward Euler projective integration scheme are two time constants, At and T. 
Given the value of the macroscopic observable at time t, a suitable "inner" timestepper (e.g. the stochastic 
simulator) is used to compute the system evolution until time t + At. Using the values of M in the interval 
(t, t + At) we estimate the time derivative of M and use it to estimate (project) the value of the macroscopic 
observable M(t + At + T) using a Taylor expansion: 

M(t + At + T) w M(t + At) + T^- . (2.1) 



dt 



(t+At) 



Hence, we compute M(t + At + T) from the value of M(t) by running the inner integrator for time At only. 
Other, more sophisticated projective integration schemes can be readily constructed [rH fisl l24| . 
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Fig. 1. Schematic of a coarse time stepper. 

If the evolution equation for M(t) is explicitly available, it is straightforward to compute M(t + At) from 
M(t) using a suitable discretization of this available evolution equation. However, if the only information 
for the time evolution of the system comes from an individual-based, stochastic model, then we have to use 
the idea of the coarse timestepper j^S] as illustrated schematically in Fig.l. Given a macroscopic variable 
M(t), we construct consistent microscopic initial conditions; we call this the lifting procedure. Next, we 
evolve the system using the microscopic simulator (e.g. kinetic Monte Carlo) for time At. Now we compute 
M (t + a, Ai) from the microscopic data for various instances < a\ < • • • < cxk — 1 (the restriction step). 
Having computed M at these instances, we use them to estimate the time derivative of M and use (2.1) as 
in the deterministic case. The resulting method is coarse projective integration. 



Some computational gain from projective or coarse projective integration can be expected provided we can 
choose At <T. A relatively large extrapolation step may save substantial computational time, considering 
the computational demands of a particle-level simulator. On the other hand, large steps can lead to low 
accuracy in simulating the dynamics of the system, or even cause numerical instabilities. We will discuss 
how one can obtain increased accuracy in projectively integrating systems with solutions evolving along 
continuous symmetry groups. The key idea is to evolve the solution (macroscopic observable) in a coordinate 
frame which tracks the evolution across space (for problems with traveling solutions) and across scales (in 
problems with self-similar solutions). 



2.1. Projective integration in a co- evolving frame 

In many cases of interest, the long time macroscopic dynamics do not involve stationary solutions but 
rather traveling, rotating, or scale invariant (e.g. self-similar) solutions I4|. Accuracy concerns in the direct 
application of projective integration to problems with such solutions [a, l3l| limit the projective time step 
T. Consider a traveling wave solution for a problem with translational invariance: it is natural to study 
its evolution in a co-traveling frame, where the solution asymptotically appears stationary (the traveling 
has been factored out). In the same sense, it is natural to study self-similar solutions in a dynamically 
renormalized frame, where the scale evolution of the solution has been factored out. Recently, a template 
based approach has been developed for the investigation of problems with translational invariance 

H3 (see 

also 0) and has been extended to the study of self-similar solutions If the description of the 

macroscopic dynamics involves scale invariant partial differential equations (PDEs), template conditions 
can be applied to derive equations describing the evolution in a dynamically renormalized framework. The 
steady state of the renormalized equations correspond to self-similar solutions of the original problem and 
the similarity exponents can also be conveniently computed |3|. This dynamic renormalization concept can 
also be applied to multiscale system models where an explicit formulation for the macroscopic evolution 
equation is not available |ct M. H2l l22l l35| . 

In this paper, our goal is to study coarse projective integration for such multiscale atomistic and/or 
stochastic problem models. To explain the idea of the co-evolving frame for such problems, it is easier to 
start with a deterministic example. We consider the PDE written in the following form 

^=C X (M). (2.2) 

Here, M = M(x,t) E B and C x : B — > B where B is a suitable Banach space of functions mapping f to E 
and the subscript x denotes the independent space variable. We define the shift operator Sc ■ B — * B and 
the rescaling operator TZa,b : B — > B by 

Sc(f) :x^f(x + C) and H A , B (f) : x -> Bf (2.3) 
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Fig. 2. (Schematic) Direct projective integration fails to produce the correct traveling shape at time t + At + T. Projective 
integration in a co-cvolving frame gives results with higher accuracy for the same time step. 

for any A, B > and Cel. We distinguish two cases - projective integration in a co-traveling frame in 
Section 2.2 and projective integration in a frame which scales with the solution in Section 2.3. 

The appropriateness of a co-evolving frame for projective integration is schematically illustrated in Fig. 2 
for a constant shape traveling wave. Direct projective integration uses the computed wave at different time 
instances to estimate its time derivative. For the instances shown in the Figure, projection according to 
(2.1) produces manifestly wrong results for large T; On the other hand, projection with the same data in 
a co-evolving frame gives results with much higher accuracy for the same time step T (i.e. for the same 
computational cost). This is because the time derivative is much smaller (here practically zero) in the 
co-evolving frame. 



2.2. Systems with translational invariance 



Let the differential operator C x in (2.2) satisfy the translational invariance property, i.e. the following 
relation holds for every Cel: 

C X S C = S C C X . (2.4) 
Let M(x,t) € B be the solution of (2.2) and let C(t) be a differentiable function of time. We define 

M(x, t) = S C ( t )M(x, t), which means that M(x, t) = S_ C ( t )M{x, t). (2.5) 
Using (2.5) and (2.2), we obtain 

dM „ — dCdM , . 

■W = c - M+ *lto- (2 - 6) 

If C(t) is given, then solving (2.6) provides the same information as (2.2). We have the freedom to choose 
C(t); we will do it so as to naturally take into account the "traveling component" of the solution. To find 



an appropriate shift C (£) , we impose an additional algebraic constraint (template condition) [^, l3^ | . The 
purpose of this template condition is to determine a shift, C(£), such that Sc(t)M is as "independent oft as 
possible." If M(x,t) were a constant shape traveling wave, then C(£) would be the change in wave position 
with time and Sc(t)M would be stationary. Hence we need a way to measure how far the wave has moved 
so we can shift it back by that amount. One can construct suitable templates in many ways. A seemingly 
natural way is to ask for the shift that minimizes some norm of the difference between the shifted wave and 
some fixed waveform (the "template", T(x) ): 

^\\S c M(x)-T(x)\\=0. (2.7) 

That fixed waveform could be M(0, x), something believed to approximate the final solution, or, in principle, 
anything else. An alternative is to use the centroid of the absolute value of the wave or some other charac- 
teristic that identifies "where the wave is" and apply a shift to bring this feature to a constant position in 
space. For a single-humped wave, one might consider using the location of the wave maximum. However, if 
during a transient the wave develops a second maximum this would clearly fail. The centroid of the absolute 
value is unique and easy to compute. If the wave is positive (or of constant sign), such as a density measure, 
the centroid has the advantage of being a linear of the wave shape that will be preserved under projective 
integration. We will formally write the template condition as the algebraic equation: 

h(M,T) = h(S c{ t)M,f) = 0, (2.8) 

where H is a functional mapping 1 x B to M. This, together with (2.6) describe the dynamics of the shifted 
solution M{x,t) as well as the dynamics of its shift, dC(£)/d£. Such template conditions arise naturally in 
the computation of limit cycle solutions in autonomous dynamical systems, where they are often also called 
"pinning" conditions |5llll|. 

If the operator C x is available explicitly, we can use (2.6) to estimate the time derivative of M at a given 
time £ and (2.1) can be used to make the extrapolation in time for simple projective forward Eulcr. To 
complete the projective algorithm in the co-evolving frame, we also must specify how the shift, C(£), evolves 
during a projective time step. As in (2.1) we approximate the shift evolution by 



dC 

C(t + At + T) w C(t + At) + T— 

dt 



(2.9) 

(t+At) 

Finally, the unshifted projected solution M(t + At + T) is computed, if required, from 

M = S_ c{t+At+T) M. (2.10) 

2.3. Systems with scale invariance 

Consider a scale invariant problem where the differential operator C satisfies the property 

CRa.b = A a B b - 1 n A , B C, (2.11) 

i.e. there exist constants a and b such that the above relation holds for every A, B > 0. Equivalently 

C x (BM Q)) = B b A a £ y (M(y)) where y = j (2.12) 

where C x and C y denotes the action of the operator C on the coordinates x and y respectively. Note that the 
system must also satisfy the translational invariance property (2.4); however, for simplicity, we will assume 
that we do not have traveling solutions, but concentrate on self-similar ones. The combination is considered 
in Section 5. 

We study solutions M(x,t) of (2.2); choosing scaling factors A (for space) and B (for the solution ampli- 
tude) as well as a reparametrization of time t(£) leads to the study of the rescaled solutions M(y,r) (see 
i) 



M(x,t) = B(r)M[j^,T(t)). (2.13) 
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Equation (2.2) can be re-writtcn in such a dynamically renormalized form as follows |3l l31| 
{ IdB-zy IdA dM dM\ dr , .\ 

[b d7 M ' Ad^ y ~dy~ + T^Jdl = AB ^ ■ (2 - U) 

Motivated by the search for self-similar solutions we select the time reparametrization r(t) as 

^ = A a B b - 1 (2.15) 

dt v ; 

which leads to 

_ = ---A/ +^V W (2-16) 

For self-similar solutions as well as -g^r are constants whose particular values depend on M. In our 

renormalization algorithm M is selected by our choice of template condition(s). 

The main idea of projective integration in a co-evolving frame is to factor out the scale evolution, so as to 
obtain a (rescaled) solution that evolves more slowly. As in the traveling case, we exploit the template-based 
approach in order to compute solutions which are "as scale invariant as possible" . A schematic description 
of the projective integration algorithm to scale invariant problems, is shown in Fig. 3. 

To determine the evolution of the scale parameters A and B we need to apply template conditions that 
control the spatial extent of the solution (A) as well as its amplitude (B). We could, for example, minimize 
the distance between the rescaled current solution and a template function to determine both A and B. 
Alternatively we could determine the amplitude by maintaining the constancy of some norm of the solution 
- the Z^-norm for a positive function would simply maintain the total mass - while for the spatial extent we 
could keep a moment of a positive measure of the solution, such as \M(x,t)\x 2 dx, constant (assuming 
that this integral is well defined) . In general, we need two independent conditions to determine the two scale 
parameters; these will take the form of two algebraic equations which are used along with (2.16) pi. lil l3l| 
to evolve the dynamically renormalized problem. These two algebraic equations take the form 

h A (M(y),T 1 ) = h A (^M(Ay),f^j =0 (2.17) 

h B {M{y),T 2 )=Q h B (^M{Ay),T^ =0 (2.18) 

where T\, T 2 are template functions. While there is considerable freedom in the choice of these algebraic 
conditions, it is important that they yield a unique and computationally simple solution for the scaling 
factors A and B. 

We can apply projective integration to the rescaled solution M in the original time frame t, or in the 
rescaled time frame r chosen above (or, for that matter, in any other convenient time variable). Using r, 
projective forward Euler is: 



— — - dM 

M(T pro j ec t) = M(t 2 ) + {Tp ro j ec t — T V~Q^ 



« M(r 2 ) + ( Tproject - r 2 ) M(T2) M(Tl) (2.19) 

T 2 - Tl 

T2 



where the rescaled times t\,t 2 and T pro j ect correspond to times t, t + At and t + At + T respectively. In 
order to approximate numerically the right hand side of (2.19) we must determine the relation between time 
t and the rescaled time r. 

We will assume that during the projection step the parameters 

k-IdP Cs = sdT (2 - 20) 

remain constant (this is true for self-similar solutions, and is analogous to assuming that the velocity is 
constant during a projection step for traveling problems). The evolution of the scale factors A,B for a 
self-similar problem is of the form 0: 

A(t)~\t-f\i, B{t)~\t-t*\ s (2.21) 
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Fig. 3. Schematic description of template based time projection for scale invariant problems. 

where t* is an appropriate (positive or negative) blow-up time, and 7, 5 are the similarity exponents. A 
typical example is the ID Barenblatt solution, the self-similar solution to the porous medium equation 
u t = (u 2 ) xx In this case, the similarity exponents are 7 = 1/3 and S = —1/3. If we consider the 

case, where the blow-up time t* — (the initial datum is a Dirac mass at the origin), then the scale factor 
A can evolve as depicted in Fig. 4(a). In Fig. 4(b), one can see the linear evolution of log A with respect to 
the rescaled time r. It can be shown [3j that the relation between rescaled time r and time t has the form 
t ~ logt; \og(B) behaves similarly. We thus expect better accuracy when the projective scheme is based on 
exponential growth of the scale factors A and B. 

During the projective step the evolution of A and B is described by: 

A(t) = A( n ) exp(U(r - n)) (2.22) 



B(t) = B( Tl ) exp(Cij(T - ri)). (2.23) 
The values A(t) = A(n),B(t) = B(ri),A(t + At) = A(T 2 ),B(t + At) = B(t 2 ) obtained through the 
application of the template conditions are used in the relation 

We can now derive an expression between the rescaled time step t 2 — t\ and the time step At. Namely from 
(2.15) and (2.22) - (2.23), we get 
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Fig. 4. (a) Evolution of scale factor A as a function of time t for the self-similar solution of ID porous medium equation. The 
similarity exponent 7 (see text for details) is equal to 7 = 1/3. The blow-up time is t* = 0. (b) Linear evolution of log A with 
respect to rescaled time r (r ~ logt). 

A(T 1 )- Q B(r 1 ) 1 - b 



-{ cxp [ - <4(r 2 - Tl ) + (1 - b)i B {T 2 - n)] - 1} 



- 1 \ = At. 



-aU + (1 - b)^ B 

It is straightforward to evaluate the parameters £4 and £b from (2.24) - (2.25) to obtain 



-aA* + (1 - b)B 
Then we compute the rescaled projection time T pro j ect 

A( n )- a B( n y- b 



{ exp [ - aA* + (1 - b)B*\ - l}(r 2 - n) = At. 



from 



{ exp [(-a^ + (1 - ft^flJCTproj-ecf - Tl)] - lj = T + At. 



-aU + (1 - &)e B 

Finally, we can also obtain the projections of the scale factors A, B from (2. 22), (2. 23): 



.4 



project 



A{T pro j ect ) — A(Ti)exp [tlA(l~project — T\) 



Bproject — B(Tproject) — B(ri)exp \_£,B{Tproject — Tl)] 

and, if desirable, recover the projection of full solution M(t + At + T) from the rescaling relation: 

x 



M(t + At + T) = B. 



project 



M 



A 



project 



' j ^project 



(2.26) 
(2.27) 

(2.28) 

(2.29) 
(2.30) 

(2.31) 



3. Systems with translational invariance - A reaction-diffusion problem 



In this section we demonstrate the efficiency of the proposed projective integration scheme in a co-evolving 
frame for a reaction-diffusio n sy stem with translational invariance. Our stochastic model is motivated by 



the Nagumo equation [jl^Eigf. 



du d 2 u 

m =D d^ +u{l - U){u ~ a) 



(3.1) 



where u denotes the reactant's concentration, a is a kinetic parameter and D is the diffusion coefficient. 
We use D = 1 in what follows. We consider a large (effectively infinite) domain with zero flux (Neummann) 
boundary conditions. The Nagumo equation is a a well known example of a parabolic system that exhibits 
traveling waves and has an explicit wave solution u(x, t) = u(x — c) given by: 

-l 



u(x) 



1 + exp 



V2J 



dc 
dt 



-V-2 (\-a 



(3.2) 
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Motivated by the Nagumo kinetics we construct a particle-based simulator, where a set of chemical reaction 
steps as well as diffusion steps are incorporated. We consider the following set of reactions: 



2N + H ^ 3A (3 .3) 

N 0. (3.4) 

The reaction rate constant for the production of reactant N is k\ = 1 + a, while the consumption rate 
constants are respectively fc_i = 1 and &2 = a. The concentration of reactant H is assumed to remain 
essentially constant and equal to 1 (H = 1). 

A standard way to simulate a spatially homogeneous chemical system is the Gillespie SSA [lj^. At each 
time step of the algorithm a pair of random numbers is generated in order to answer two essential questions: 
when will the next event - chemical reaction - occur and which reaction will it be? We incorporate in our 
system the effect of spatial diffusion too; N diffuses with a diffusion coefficient, D. 

The generalisation of Gillespie ideas to spatially distributed systems can be found in e.g. Here, 
diffusion is treated as another set of "reaction steps" in the system. The domain of interest is discretized into 
J lattice sites with constant distance h between them. We denote by Ni the number of respective molecules 
at lattice site i. This means that we describe the state of the stochastic reaction-diffusion system by a 
J-dimensional vector N = \Ni , N2 ■ ■ ■ , Nj] , and the following reactions at each time step are considered 

2N t + H ^ 3Ni 

*-i } i = l,...,J, (3.5) 

N J% 

Ni -^N i+1 , z = l,...,J-l, (3.6) 

Ni -1+Ni-!, i = 2,..., J. (3.7) 
The set of reactions (3.5) implies that the reaction mechanism (3.3) - (3.4) is implemented at each lattice 
site of the domain. Moreover, diffusion is introduced as a set of new reactions (3.6) - (3.7), whose transition 
rates are denoted by d. The transition rates for d are connected to the macroscopic diffusion coefficient, D, 
which at a certain limit (h <C 1), is given by the formula d — D/h 2 . The augmented set of reactions (3.5) - 
(3.7), together with suitable boundary conditions, can thus be simulated using Gillespie SSA. 

We will start by projectively integrating the Nagumo partial differential equation, and then we will 
illustrate the coarse variant of the method for the particle-based implementation of the scheme (Gillespie 
algorithm) . 



3.1. Nagumo Equation - PDE description 



We first consider the deterministic description (PDE) of the Nagumo problem and illustrate the accuracy 
improvement to the projective method from operating in a co-evolving frame. In our numerical computa- 
tions a = 0.01, so that velocity of the traveling wave is dc/di w —0.693 according to (3.2). The (long) 
one-dimensional domain [—30,30] is discretized into 601 equidistant nodes (i.e., the distance between two 
successive nodes is Sx — 0.1). The spatial partial derivatives are approximated with central finite differences 
and the applied boundary conditions are of Ncumnann type. The initial condition is: 

-30 < x < 0, 

0<x<10, (3.8) 
10 < x < 30. 

The time step of the inner integrator (here a simple forward Euler) is 6t = 10~ 4 to satisfy the stability 
criterion 28t < Sx 2 . 



j for 

u(x,0) = < x/10 for 
I 1 for 
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Fig. 5. Nagumo example: Solution obtained from direct simulation and (a) non co-traveling, (b) co-traveling projective inte- 
gration at time t = 15 (see text). 

A typical projective integration step requires the solutions u\,u 2 at two distinct reporting times ti,t 2 . 
To obtain the solution at projection time t pro j ec t, we simply apply the Taylor expansion (2.1). We choose 
two reporting times t\ = 0.1 and t 2 = 0.2 which correspond to 2 x 10 3 steps of the inner integrator. A 
projection time t pro j ect — 0.5 thus saves 3 x 10 3 inner integration steps. If we use a projective method 
which ignores translational symmetry the results are manifestly inaccurate. Fig. 5(a) compares the results 
of projective integration to those of full direct simulation; taking translational invariance into account (see 
Fig. 5(b)) clearly shows the improved accuracy. 

The template condition used to obtain these results was: 



/30+c 
-30+c 



u{y)dy 



/30 <-30 
u(x + c, t)dx = / u(x,0)dx, 
-30 J-30 



(3.9) 



implying that the integral of the shifted solution remains constant and equal to the integral of the initial 
condition in the domain of interest. Application of (3.9) at each reporting time t\, t 2 produces the "shifted" 
solutions iii, U2 and the corresponding shifts ci, c 2 . The projection of u is obtained from 



^{tproject) — 11*2 T" [tproject ^2 ) ~7T7 



dt 



u 2 + (t 



project 



h) 



u 2 - u 1 

t 2 -h 



and the projection of the shift c from 



^■{tproject) c 2 ~t~ [tproject t 2 ) 



C 2 + (t 



project 



-*2) 



1-2 



c 2 - Ci 
ti-tx 



(3.10) 



(3.11) 



Finally, the full projected solution, reconstructed in physical space x is recovered, if desired, applying the 
inverse shift operator: 

u(x, tproject) — u(x c{t pr oject) i tproject) ■ (3.12) 

The benefit of projecting in a co-traveling frame is more clearly depicted in Fig. 6, where we plot the error 
evolution (L 2 norm of the difference between the solution obtained from direct simulation, Udirect and the 
solution computed from projective integration (co-traveling or not), upj at the same time, t), i.e.: 

1/2 



e(t) 



V-30 



(Udirect{x,t) - U PI (x,t)) 2 dx 



(3.13) 



The error e(t) resulting from the non co-traveling projective method increases with time, while the results 
of the application of the modified projective integration scheme appear highly accurate. The increased 
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Fig. 6. Error evolution for the non co-traveling and for the modified co-traveling projective integration in Fig. 5. 




time t 



Fig. 7. Time derivative of the computed Nagumo solution (a) u and (b) u up to t = 15. The relatively high value of ut is the 
primary cause for low accuracy of the obtained results when projective integration is applied in a stationary frame. 

accuracy of the modified method is due to the slow evolution (compared to the faster evolution of the 
unshifted solution) (see Fig. 7). 

For the co-traveling computations the solution after some time appears stationary (Fig. 8(a)); this is (an 
approximation of) the stable Nagumo traveling wave. Its (constant) speed (Fig. 8(b)) can be approximated 

by 

J « (3.14) 
dt t 2 -ti 

3.2. Coarse projective integration in a co-traveling frame - Kinetic Monte Carlo simulation of a 
reaction- diffusion system 



We now apply the same methodology to traveling problems for which the model simulations are conducted 
at a microscopic (stochastic, particle) level. 
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In our illustrative example the particle-based simulation is the Gillespie SSA presented in Section 3 applied 
to the same kinetic scheme. The one-dimensional domain of interest [—30, 30] is discretized with J = 601 
lattice sites. The distance, h, between two successive lattice sites is equal to, h = 60/600 = 0.1. The zero flux 
boundary conditions are incorporated applying a zero reaction rate for the "reactions" (3.6) and (3.7) at sites 
i = J and i = 1 respectively, i.e. we do not allow the particles at i = J to diffuse to the right and particles 
at i = 1 to diffuse to the left. At the deterministic limit the reaction rate constants are k\ = 1 + a = 1.01, 
k-i — 1 and &2 = ct = 0.01. The macroscopic diffusion coefficient, D = 1, corresponds to a diffusion 
rate constant d = 1/h 2 . In our computations, we assume that the number of particles corresponding to 
dimcnsionless density, u = 1, is Nq = 1000. The reaction parameters for the Gillespie code have been chosen 
consistently. The number of particles at site % is denoted by N i: i = 1, . . . , 601. The initial condition is: 








for 


1 < i 


< 201, 


N t = < 


lOOi 


for 


202 < 


i < 401 




1000 


for 


402 < 


i < 601 



The results obtained from the kinetic Monte Carlo simulation are illustrated in Fig. 9, where one can clearly 
see the formation of a (stochastic) traveling interface sweeping the one-dimensional domain. The stochastic 
simulation described above can be computationally intensive, especially if one increases the number of 
particles. Such computations can be accelerated through a coarse projective integration scheme; translational 
invariance at the coarse (concentration field) level should then be taken under consideration for more accurate 
results. 

We now describe the modified coarse projective integration scheme applied to the Nagumo kinetics- 
motivated, kinetic Monte Carlo simulator. The SSA is performed on J = 601 lattice sites. In order to obtain 
a less noisy distribution we estimate a smoothened distribution in n = 101 nodes, using the local averaging 
operator 



1 fe=4 1 fe=6i-3 , fe=601 

Mi = -^^, M i = - Yl N k, i = 2,--- ,100, M wl = -Y, N k- (3-16) 

k=l fe=6i-9 fe=598 

We would like to approximate this distribution in Fourier form; due to the boundary conditions, we 
consider the difference distribution / (in effect, the spatial derivative) defined as: 
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Fig. 9. Gillcspic-bascd time evolution of a reaction-diffusion system motivated by Nagumo kinetics 



and its Fourier approximation [27 



fj = 
, i.e.: 

2 



Ma 



Mj-0 for j = l, 
M,-_i for j = 2, ..,n 



fe=K 

E 

fc=i 



at- cos 



27TX 



bu sin 



27TX 



(3.17) 



(3.18) 



where L is the domain length (L = 60). The coarse variables in our computations are the first K Fourier 
coefficients of /. In the projective integration context we compute these Fourier coefficients at two reporting 
times t\,t2, approximate their time derivatives at tj, and extrapolate to the projection time t pro j ect . Such a 
computation, does not take into account the translationally invariant character of the problem, leading to low 
accuracy results for relatively large steps. Application of the coarse projective scheme in a co-traveling frame 
can capture the dynamics of the same system with enhanced accuracy, even for relatively large projecting 
horizons, T = t project - i 2 . 

We denote the shifted version of /, with /, i.e.: 

/(*) = /> - c). 

The Fourier coefficients &j, bi of / are then given by: 



a, cos 



We apply the template condition 

r 30 



\2ttc~ 




' 2?rc" 




+ bi sin 









" 2ttc" 




' 2?rc" 




+ bi cos 




_ 







d r u - ~ 



Ida; = 



d r 30 

— / f(x + c)T(x) dx = 
dc /-30 



(3.19) 
(3.20) 

(3.21) 



seeking maximum overlap between the shifted solution / and a template function T. Choosing the trigono- 
metric template function T(x) = 1 — cos[27ra;/L] reduces the template condition to: 



arctan 



(3.22) 



dc L 

The computation of the shifts c\,C2 at reporting times ti,i2j enables the determination of the "shifted" 
Fourier coefficients Si, bi and their projection at time t pro j ec t'. 

&i(t 2 ) - &i(ti) 



Cli(t 



project } 



project 



t 2 )- 



for 



0,...,K 



(3.23) 
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^{(tproject) — ^(^2) (tproject ^2) 



6«(*2) - 



for i = 1, X. 



(3.24) 



*2-*l 

The projected / is then recovered by applying (3.18). The distribution M of particles at the n nodes is 
then computed and the lifting procedure concludes by interpolating M to the J lattice sites of the kinetic 
MonteCarlo time simulator through MATLAB's intrinsic function interpft. When this interpolation does 
not give an integer number of particles we round off. This gives us the projected particle distribution in 
the co-evolving frame. The spatial position of this distribution is determined by the projection of the shift, 
Cproject, according to 



Cproject — C2 4" (tproject ^2) 



C2 - Cl 



(3.25) 



*2-*l 

Note that the projected solution does not necessarily exactly satisfy the template condition; we circumvent 
this issue by calculating the 61 (tproject) from (3.22) with c — Cproject and (i\ — d\ (tproject) • 
The Fourier coefficients of the full, un-shifted solution, / are 



di(tproject) — (tproject) COS 



. 2ttc 



■project 



bi(t 



project ) S1H 



. 2ttc 



project 



(3.26) 



bi{tproject ) Oj (t 



project 



) sin 



. 2-kc. 



■project 



L 



+ h{t 



project 



) COS 



. 2ttc. 



project 



L 



(3.27) 



from which the projected particle distribution at the j lattice sites can be obtained. 

Results of this template- based projective scheme are presented in Fig. 10, accurately capturing the (coarse) 
dynamics of the kinetic Monte Carlo Nagumo simulator, even for relatively large projection steps. The 
reporting times at each projective step were taken so that t% — t\ = 0.25, the projection horizon was 
tproject — t2 — 0.5 and the number of Fourier coefficients was K = 15. The coarse variables of the problem 
(Fourier coefficients) become essentially constant in the co-traveling frame (see Fig. 11). 

4. Projective and Coarse projective integration in a dynamically rescaled frame: Diffusion 



In this section we study a projective scheme modified for scale invariant systems, in particular systems 
that possess self-similar solutions. Our illustrative example is simple one-dimensional diffusion, both as a 
deterministic PDE and via a Monte Carlo-based simulation. 
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Fig. 11. Nagumo problem. Evolution of coarse variables ((a) 3rd Fourier coefficient 03 and (b) 4th Fourier coefficient 124) in the 
co-traveling frame (stars) and in a constant frame (dots). 

4.1. One-dimensional diffusion - PDE example 
We study the simple mass diffusion equation 



u t = u x 



(4.1) 



It possesses well-known self-similar solutions; we will exploit this property, in order to perform relatively 
large projective steps accurately. One can easily verify the scale invariant character of (4.1) 



C x (Bu (£)) = BA- 2 C y (u(y)) where y = £ 



(4.2) 



For our numerical computations we discretize the one-dimensional domain x € [—10, 10] in 1001 equidistant 
nodes (i.e., dx — 0.02). The dynamically renormalized diffusion equation (along the lines of (2.16)) is 



du „ 1 dB „ 1 dA dii 



(4.3) 



The spatial derivatives are approximated with central finite differences and we consider zero flux boundary 
conditions. The selected initial condition is 



u(x, 0) = 



for \x\ > 1, 

1 for Id < 1. 



(4.4) 



At each step of the projective integration scheme, we choose two reporting times ti, t<i and the solutions 
there, u\, U2 respectively. 

A co-evolving frame formulation requires rescaling of the computed solutions using 



u(x, t) — B(t)u 



Mr) 



(4.5) 



as discussed in Section 2.3. We can evaluate both the rescaled solution, u, and the scale factors A, B at each 
reporting time step by solving two template conditions. In this illustrative example the template conditions 
chosen are 



/+oc r+00 
u(y)T 1 (y)dy = / -u^T^dy = 



(4.6) 
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Fig. 12. Projective integration in a co-evolving (co-collapsing) frame applied to the one-dimensional diffusion PDE example, (a) 
Instances of the evolution of rescaled solution u obtained at different -rescaled- projection times r vro j £ct . The rescaled solution 
u converges to a steady state profile, corresponding to a member of the self-similar family solutions. The initial condition is also 
depicted in the figure, (b) Scale parameter £a = d 1 °^ A ^ and £s = d lo s( B ) 
to the procedure described in Section 2.3. 



values computed at each projective step according 



for the template function 



Ti(y) = 



-1 for \y\>l/2, 
1 for \y\ < 1/2 



and 



u(y)T 2 (y)dy = n 



—u(Ay)T 2 (y)dy = fi, 



(4.7) 



(4.8) 



where fi is a constant and the second template function is chosen as T 2 {y) = 1. The template condition 
(4.8) keeps the mass of the rescaled system constant, equal to the initial mass. Below we present results in 
a co-evolving frame, for t 2 — h = 0.1 and a projection step of t pro j ect — t 2 = 0.2 (thus economizing 10000 
time steps of the inner Euler integrator). The evolution of u is depicted in Fig. 12; a stationary profile is 
approached after some "rescaled" time, r; this profile is a member of the family of self-similar solutions of 
the diffusion equation. The scale parameters £4 = d 1 °^ A ^ and £g = HsiLjE) shown in the same figure also 
approach stationarity. 
Comparison of the errors 



e(i) 



10 



("direct fa, t) - Upi{x, t)f Ax 



(4.9) 
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for projective integration in a co-evolving frame with those for unmodified projective integration (shown in 
Fig. 13) illustrates the advantage of projecting in a dynamically renormalized frame; Udirect is the solution 
obtained from direct simulation of (4.1). The errors are computed for the reconstructed solutions upi. 

Once more, the extra accuracy can be attributed to the slower evolution in the dynamically renormalized 
frame (see Fig. 14). 



4.2. Random walker simulation of one- dimensional diffusion - Coarse Projective Integration 



In this section we present "renormalized coarse projective integration" applied to a Monte Carlo algorithm 
simulating diffusion in a population of 10 6 random walkers in one space dimension. 

The macroscopic observable in this case is the cumulative distribution function (CDF) of the particle 
positions denoted by /. The domain of interest [—10, 10] is discretized into 1001 equally spaced nodes 
(dx — 0.02). The CDF is then determined as a function of the discretized spatial domain {x\, ziooi}- 
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Fig. 13. Comparison of error evolution for the simulations in Fig. 12 when the projective integration and the modified projective 
integration is applied. The modified projective integration algorithm application manages to produce accurate results even at 
the early stages, where the solution still evolves towards its self-similar shape. 
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Fig. 14. One-dimensional diffusion equation, (a) "Time" derivative of rescaled solution u at T2-reporting times, using the (4.6) 
and (4.8) template conditions, (b) Time derivative of u as obtained from the application of the original projective integration 
scheme at t2-reporting times. The relatively high values of time derivative ut is the main reason for the failure of the method 
to capture the correct dynamics of (4.1). 

Each particle, i, is described by its position X,. The Monte Carlo time step is St — 0.0001 and during each 
time step each particle will randomly move left or right with equal probability by an increment 5X = \/2St. 
Denoting CDF at mesh point x t as fa , the probability density function of particles is evaluated by 

fi Ji— 1 



N, 



1/2 



(4.10) 



where -/Vj_i/ 2 is the macroscopic density of particles at the midpoint [xi-\ +Xi]/2. Lifting - i.e. construction 
of a microscopic state consistent with density (4.10) - is done as follows. The particles are placed in space 
so that their density piecewise linearly interpolates the midpoint values (see Fig. 15). 
We assume that an evolution equation for the CDF of particle positions f(x) exists: 

dt 



£*(/)■ 



(4.11) 
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Fig. 15. The piecewise linear distribution of density JV corresponds to a microscopic realization of particle positions. The 
distribution N is evaluated as the spatial derivative of the CDF. 

Before applying the coarse projective scheme in a co-evolving framework we should test the scale invariance 
of the unknown differential operator C x and extract its scaling exponents. In our computations there is no 
amplitude scaling since / is a CDF. The operator should satisfy 

C x (f(^j)=A a C y {f{y)) where y = j (4.12) 

for any A. Despite the fact that the explicit formulation of C x is unknown we can estimate its action on test 
distribution f(x) through the computation of We perform short computational experiments to estimate 
the action of operator C x as follows: 



(1) We select a test function ipo and a positive constant A. 

(2) We initialize the kinetic Monte Carlo simulator so that the CDF of particles is equal to po (using 
the lifting procedure in Fig. 15). We run the kinetic Monte Carlo simulator for a relatively short time 
interval (in macroscopic terms) ST. We obtain the new CDF ipi and estimate the time derivative of 
ip from the expression 

^ « (4.13) 
dt ST v ' 

(3) We initialize the kinetic Monte Carlo simulator so that the CDF of particles is equal to po(x) — 
(fo(x/A) (using the lifting procedure in Fig. 15). We run the kinetic Monte Carlo simulator for a short 
time interval ST. We obtain the new CDF (p\ and estimate the time derivative of ip similarly as in 
(4.13). 

(4) We estimate the value of exponent a in (4.12) by minimizing the residual 



R{a) 



dip 




A dt 




dt 


Ax 


X 



(4.14) 



where 



denotes the standard Euclidean norm. 



In this case we use: 

where po^ is the value of test function tpo at mesh point Xi, B(7, S) is the Beta function with parameters 
7 = 8 and 5 — 10, ST = 0.01 and A = 1.15. The results are shown in Fig. 16. The procedure described 
above can be performed for different values of A and different test functions po . The value of the exponent 
a minimizing the residual (4.14) was in all tested cases close to —2. It confirms the operator's C x scale 



invariance property ( 



Ax 



almost coincides with A a ^ 



). The value of the scaling exponent a = —2 is 



used in the proposed co-evolving projective integration scheme. 
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(a) 



(b) 



Fig. 16. (a) Time derivative of test function ipo given by (4.15) (solid line) and the rescaled test function ipo(x) = ¥>o(-j) 
(dashed line); parameters of algorithm (1) - (4) arc given in the text, (b) Testing the scale invariance property of operator C x 
(see (4.12)) scale invariance property. The value of scale exponent is a = —1.97 ~ —2. 



In our modified coarse projective scheme we evaluate the macroscopic observables at k = 2 distinct 
reporting times ti, t<i with corresponding CDFs /i, fi- At each step of the projective integration scheme the 
time step is — t\ = 0.05; the projection step is t pro j ec t — ti = 0.1. The initial CDF fo at mesh point Xi is 
given by 



h{xi) = fo,i 



for 1 < i < 451, 
451]/100 for 452 < i < 551, 

1 for 552 < i < 1001. 



(4.16) 



Both the scale factor A and the rescaled solution /, where 

f{x) = f(Ax), 



(4.17) 



are obtained from the application of the template condition: 



/(C2) - /(Ci) = " f(M2)-f(Mi) = ", (4-18) 

where Ci, C2 are given real numbers and v is constant. This template condition, enforces a constant number 
of particles in interval [Ci,C2]- Let us note that the CDF is defined only at mesh points {xi, a;iooi}- 
Whenever template condition (4.18) requires values of / outside mesh points {xi, ...,£1001}, we use linear 
interpolation. For our computations we chose £1 = —0.25 and £2 = 0.5, while the constant v is evaluated 
from the initial condition fo, i.e. /o(C2) — /o(Ci) = v - 

As in the PDE diffusion example, we choose to evolve the rescaled CDF in rescaled space y and time r. 
The temporal Taylor expansion is performed in terms of r; we therefore need to evaluate the n, t 2 and 
Tproject values corresponding to t\, t 2 , t pro j ect . The procedure was reported in Section 2.3; note that here 
the B scaling is omitted since / is a CDF. 

The numerical results are shown in Figures 17 and 18. In Fig. 17 we present results of the first five coarse 
projective steps applied to the rescaled CDF /. The circle-marked lines correspond to reporting ri-times, 
the square-marked lines to reporting T2-times and the triangle-marked lines to projective T pro j ec t -times. In 
Fig. 18 we plot the evolution of scale factor A both in terms of time t and of rescaled time r, as computed 
from the modified projective integration using the template condition (4.18). 
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(a) 



(b) 



Fig. 17. Random walker simulation of one-dimensional diffusion: (a) Coarse projective integration applied to the CDF / in the 
dynamically co-evolving frame. The lines marked with circles correspond to solutions, the square-marked lines correspond 

to /(T2) solutions and the lines marked with triangles correspond to / obtained from projection at T pro j ect . (b) Time derivative 
of / evaluated at the I s '— 5" 1 projective steps of the co-evolving projective integration algorithm. 
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Fig. 18. Random walker simulation: Scale factor A values computed from template condition (4.18) during the application of 
the modified projective integration algorithm. The circle points correspond to A-values computed at (a) t\ - reporting times 
(b) ti - reporting times. The square points correspond to (a) t2 and (b) T2 = i~(t2) reporting times. The A values obtained at 
projection times t pro j ect ((b) Tp ro j ect ) are marked with triangles. 

5. A general approach for problems with asymptotic or approximate scale and translational 
invariance 



In most real- world problems the issue of translation and rescaling must be handled simultaneously. Even 
if the problem is self-similar, our choice of templates for normalizing may lead to an apparent translation 
with time. For example, the viscous Burgers equation 

U t = UU X + KU XX (5.1) 

has a self-similar solution 

u(x,t) =t-^ 2 w(xt- 1/2 ) (5.2) 

where w(y) is not a symmetric function. Unfortunately if we choose the wrong origin in the y coordinate, 
we will find that the evolving waveform is also traveling because we will actually be looking at 

u(x, t) = t-^wdx - x )t- 1/2 + xot- 1 / 2 ) (5.3) 
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where the x t 1 / 2 looks like a translation with time. 

Many problems are only asymptotically self-similar. For example the solution of 

u t = (1 + u 2 )u xx (5.4) 

asymptotically approaches the self-similar solution of the heat equation because as u decays the u 2 term 
becomes asymptotically small. Other problems may be approximately self-similar. For example, they may 
approximately satisfy a scaling relationship such as 

Cx ( Bf (j)) = B b A°C y (f(y))[l + 0( £ (t))] (5.5) 

where e(t) is small and y = x/A. 

Therefore, we want to look for time-dependent translations and rescalings that yield a slowly varying 
waveform that can be integrated in time more accurately because of its smaller time derivatives. However, 
we cannot use the mechanism described in the previous sections, where we examined "purely" scale invariant 
systems, because we do not know a and b and they may not be defined away from the asymptotic limit. 

As the integration proceeds, it generates successive values of the solution, u(x,t). Periodically we need to 
apply three transformations to u(x, t) to get 

u(y,t) = ^--u(C(t) + A(t)y,t) (5.6) 

with y = [x — C(t)]/A(t) to try to get a u(y 7 t) which is as independent of time t as possible. Motivated by 
the approach for co-evolving systems that are exactly self-similar, we will also consider a transformation of 
t to r(t) in the expectation that A(t) and B{t) will evolve approximately exponentially in r so that linear 
projection of their logarithms will provide an accurate approximation. 

Template conditions are needed to determine the three scalings in (5.6). If we had a moderately good 
approximation to the final waveform, it would be tempting to ask that the scalings be chosen to minimize 
the difference between u(y, t) and that approximate waveform. One might even choose the scalings so as to 
minimize some norm of the time derivative. However, such conditions are nonlinear, and nonlinear conditions 
cause a problem in the projective step. A projective step takes the form 

, tp ro j ec t £2 r - ~ i /r rr\ 

Uproject = U 2 H [U 2 -UlJ. (5.7) 

t2 — II 

In other words, u pro j ect is a linear combination of iii and u\. If ui and u\ satisfy a linear template condition, 
Uproject will also satisfy that condition automatically. That is not necessarily true for a non-linear condition; 
after a projection step one would have to re-apply the condition. 

Therefore, we choose linear conditions to determine the scalings. If the solution is non-negative - as it will 
be for many physically based problems in which the variables we will use in the templates are quantities such 
as density (represented at the microscopic level by numbers of particles) there is a straightforward recipe. A 
shift can be determined by demanding that the center of gravity be shifted to a fixed position, typically the 
origin. This is a trivial calculation. For example, for uniform particles we simply average their positions. In 
a continuum model we ask that 

/oo 
yu(y)dy (5.8) 
-oo 

which means that 

f _ xu(x)dx 
J_ oc u(x)dx 

The A scaling can be calculated by requiring that a certain fraction of the integral of u lies within a specified 
interval, for example that j_ x u(y)dy — 0.5 J™ u(y)dy, although it is computationally easier to specify the 
second moment and require that 

f°° y 2 ii(y)dy 

J-oo u (y) d y 
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which implies that 



.2 _ JZo x2u (x + C)dx 



J-oo u \ x )dx 

Note that the shift has been applied before the next moment is calculated. The amplitude scaling B can be 
calculated by requiring that the total mass, 



J — ( 



u(y)dy, (5.12) 

J — oo 

remain constant, say equal to fi, leading to 

f°° u(x)dx 

B = J -°° / . (5.13) 

As the microscopic integration proceeds we calculate the values of A, B and C at selected times to, t\, 
t 2 , ... and the rescaled solution u(y, to), u(y, t\), u(y, t 2 ), • • • in the coarse variables. We can decide whether 
it is appropriate to apply a projective step to the coarse solution based on the local behavior of the coarse 
variables and the scaling values. 

As we approach the region where the solution is close to self-similar, we need a way to compute the 
transformation from t to r so that we can use exponential projective integration in r. We can do this as 
follows. We assume a form like (5.5) in the PDE 

u t = C x (u) (5.14) 

and assume a solution of the form 



u(x,t)=B(r)w^—,T) (5.15) 
for some r(t) and w(y,r) which is slowly changing in r. Then we get the equation 

(lf W ~ lt yWy + WT )m= B AO/ >[l + 0(e(t))]. (5.16) 

We want to choose A(t), B(t), and r(t) (which are completely at our choice) so that if there exists an 
approximately self-similar solution, that is, a solution of 

{b l w-a 1 yw v )—=B h - 1 A a C v w (5.17) 

for some a\ and 61, w tends to it. We naturally choose r t = cB b ~ 1 A a for some constant c and A and B 
so that A T I A and B T /B are nearly constant and equal to a\ and 61 respectively. (In practice, we will be 
choosing A and B to account for the observed growth in width and amplitudes of the computed solution.) 
If A T /A and B T /B are constant, then A(t) = exp(a + air) and B(t) = exp(6 + 61 t). Hence 



dt f d T \ cx p[~ a ( a o + a i T ) — (b — l)(&o + °i T ) 
~dr~ ~ \~di 



(5.18) 



or 

expf— ana — bo(b — 1)1 . . , . ,„,„s 

t = t c - ^ a +biib °\ ))c n eM-(a.a + b l( b- 1 M (5.19) 

Note that we have this exponential behavior for t(r) regardless of the actual values of a and 6. Since the 
scale (and origin) of r are arbitrary (we are picking them), we can rewrite this equation as 

t = t c + /3exp{r). (5.20) 

Suppose now that we perform a calculation starting at t and integrate to t\ and t<i, computing A(ti) and 
B{ti) as we proceed using template conditions. We can assume that r = since the origin is arbitrary. We 
need to find n and r 2 . We have from (5.20) 

(U - to)/P = exp(r 4 ) - cxp(ro) (5.21) 
22 



or, using t = 

n = log(l + (U - *„)//?). (5.22) 
Under the assumption that A T /A is constant we have 

(and a similar relation for _B). Using (5.23), (5.22) and to = tq = 0, we get 

This can be solved for (3 and then we can use (5.22) to find the Tj. In the projective step, a suitable 
representation of the rescaled solution, w can be projected in r, and A and B are projected exponentially 



5.1. An asymptotically scale and translationally invariant PDE example. 

Below, we present some representative results from the application of this general approach to the equation: 

u t = k(1 + u 2 )u xx + uu x with k = 0.025 (5.25) 

which asymptotically approaches the solution of the viscous Burgers equation, because as u decreases the 
u 2 term becomes asymptotically small. Both translations and rescalings have to be incorporated for this 
problem as one can see from the time evolution of u in Fig. 19. Three template conditions are applied for 
the evaluation of the shift C, and the scale factors A, B at each reporting step, which have the form of 
(5.9), (5.11) and (5.13) respectively. The constants K and fi appearing there are computed from the initial 
condition, i.e.: 

f°° x 2 u(x,0)dx 

K = J ~Z / i ' (5.26) 
J_ oo u(x,0)dx 

and the initial mass: 

/•OO 

li= / u(x,0)dx. (5.27) 



The initial condition in our computations is u(x, 0) = cxp(— x 2 ). The applied boundary conditions are of 
Neummann type, the one-dimensional computational domain [—10,10] is discretized with 1001 nodes and 
the spatial derivatives of (5.25) are approximated with central finite differences. 

The direct time integration is performed with an explicit Euler scheme, with alt = 10~ 5 ensuring the 
stability of the integration for the given discretization. The three reporting times are chosen so as At = 0.1, 
while the projective step is taken to be T — 0.2. At the initial stages of the integration, the solution is far 
from its asymptotically self-similar shape solution and it is trivial to show that (5.25) is not scale invariant, 
at least for large enough values of u. In this case we can project linearly in time the shift factor C the 
scale factors A, B and the renormalized solution u which is derived from the rescaling equation u(y) = 
l/Bu(Ay + C). When the solution approaches the self-similar regime, then we can apply the transformation 
of t to t and follow the procedure described in Sec. 5. The evolution of the factors A and B computed from 
direct simulations and from the 3-step, template-based projective scheme are depicted in Fig. 20. Finally, we 
illustrate the accuracy of this method, presenting the solution computed from the direct simulation at time 
t = 9.9 and the one derived from the projective method at the same time (see Fig. 21). 

6. Summary and Conclusions 

In this paper we have illustrated projective and coarse projective integration in a co-evolving frame for 
problems with continuous symmetries, and in particular for problems with (coarse) scale invariance and 
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Fig. 19. Direct simulation of (5.25) with initial condition u(x,0) = exp(— x 2 ). The solution u{x,t) travels both across scales and 
in space, forming a steep interface which moves to the left part of the one-dimensional domain. Instances of u are presented at 
t = 0,1,3,5, 7, 10. 

(coarse) translational invariance. The system temporal evolution is observed in a traveling and/or dynami- 
cally renormalized frame, in which transient solutions approaching traveling waves or self-similar solutions 
appear slowly changing, and can be integrated more accurately because of the smaller time derivatives. 
Larger extrapolation time steps can thus be applied without degrading the accuracy of the projected so- 
lution. The simplest projective algorithm (projective forward Euler) was illustrated; more sophisticated 
multistep and even implicit projective algorithms are also possible (e.g. [Til I24I l30|V We have illustrated 
several representative examples of template (pinning) conditions that are used to dynamically define the 
coevolving frame in which projective integration takes place. Our model examples included both contin- 
uum and microscopic-based implementations. The translationally invariant, co-traveling case was illustrated 
through the Nagumo reaction-diffusion PDE |f| [2(J [2^ in one spatial dimension, as well as an SSA-based 
|l9l | stochastic implementation of the Nagumo kinetics for coarse projective integration. The scale invariant 
case was illustrated through simple one-dimensional diffusion: projective integration of the PDE version and 
coarse projective integration of a stochastic implementation involving a large ensemble of random walkers 
were presented and the results compared with direct, full simulation. Finally, we described a more general 
projective method designed for systems with asymptotic or even approximate invariance, and where scale 
invariance and translational invariance co-exist. 

The thrust of the paper was in describing and illustrating the methods, providing some evidence and 
qualitative justification for the resulting computational savings. This constitutes only the starting point for 
the numerical analysis of the algorithms, both for the deterministic and for the stochastic cases, which is the 
subject of further research, ft is worth reiterating that the template based approach transforms traveling 
or self-similar problems into steady-state ones; traveling wave speeds, "scale velocities" ^a^b and similar- 
ity exponents are simple and natural byproducts of the approach. Accelerating the computation of coarse 
self-similar shapes and coarse similarity exponents for microscopic/stochastic simulators can be useful in a 
variety of disciplines, ranging from microhydrodynamics to core collapse in star clusters |35| . In the stochas- 
tic case, the accurate and efficient estimation of time derivatives from stochastic simulations becomes a vital 
component of the algorithm, and one must move beyond simple differencing and least squares estimators 
(like the ones we used here) to maximum likelihood ones (see e.g. 0). It is worth noting that -whether in the 
deterministic or in the coarse-grained case- it is important to explore the relation between modern adaptive 
mesh techniques used for the computation of self-similar solutions [6|, (jj |29j with the template-based ones 
presented in | r J.l3l| and exploited here. The most important factor in the success of coarse-grained projective 





24 




-2.5_ L 



2.5 r 



0.5L 



4 6 
time t 



10 




direct simulation 
modified projective 



4 6 
time t 



8 



10 



direct simulation 
modified projective 




(a) 



1.6r 



(b) 



B 



0.9 
0.8 
0.7 
0.6 
0.5 



(c) 



2.5 



2.5 



direct simulation 
modified projective 





direct simulation 
modified projective 



3 

time t 



3.5 



direct simulation 
modified projective 



3 

time t 



3.5 



Fig. 20. Computed values of (a) shift factor C and scale factors (b) A, (c) B from the application of template conditions 
(5.9), (5.11) and (5.13) respectively, to the solution u(x,t) of (5.25). The solid lines correspond to results derived from direct 
simulation of (5.25). The C, A and B values, computed from the 3-step template based projective method (see Section 5.1), 
are depicted by dotted lines. 
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Fig. 21. Evolved solution u of (5.25) at time t = 9.9. The solid line represents the solution derived from direct simulation and 
the dashed line corresponds to the solution computed by the 3-step template-based projective scheme described in Section 5.1. 

integration (and of equation-free computation in general) is the lifting step: the ability to construct ensembles 
of fine-scale realizations consistent with a given macroscopic description. For the problems discussed in this 
paper, the coarse observables in terms of which the process was modeled allowed for a relatively easy and 
computationally inexpensive lifting. If the coarse-grained model is cast in terms of different coarse-grained 
observables (e.g., particle pair correlation functions) the lifting step may become much more difficult and 
expensive (see e.g. H^|). Clearly, the cost of the lifting step (a very much problem-dependent feature) must 
be factored in when evaluating the potential savings of coarse projective integration. We close by reiterat- 
ing that what we have presented is only a first step in the study of coarse-grained projective integration 
algorithms for systems with (coarse) continuous symmetries; the potential benefits illustrated here for two 
simple model problems argue that the algorithms both on the continuum front (e.g. projective integration 
for differential-algebraic equations [l^, relations to adaptive mesh algorithms) and on the stochastic front 
(issues of lifting and estimation) warrant extensive further study. 
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